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Abstract 

The precise determination of the arrival direction of cosmic rays is a fundamental prereq- 
uisite for the search for sources or the study of their anisotropies on the sky. One of the most 
important aspects to achieve an optimal measurement of these directions is to properly take 
into account the measurement uncertainties in the estimation procedure. In this article we 
present a model for the uncertainties associated with the time measurements in the Auger sur- 
face detector array We show that this model represents well the measurement uncertainties 
and therefore provides the basis for an optimal determination of the arrival direction. With this 
model and a description of the shower front geometry it is possible to estimate, on an event by 
event basis, the uncertainty associated with the determination of the arrival directions of the 
cosmic rays. 



1 Introduction 



The nature and origin of ultra high energy cosmic rays (UHECR) is a continuing puzzle of modern 
astrophysics [Tj. Essential to the study of their origin is the measurement of their arrival directions 
with a careful evaluation of their measurements uncertainties. 

The Pierre Auger Observatory ||2l is designed to measure the flux, arrival direction distribution, 
and composition of the UHECR. In order to have full sky coverage, it will have two components, 
one in each Hemisphere. The Southern site of the Observatory is in its last stages of construction 
in Malargiie, Argentina, and has been acquiring data since 2004. When completed, it will be 
composed of 1,600 water Cherenkov detectors (WCD) spread over 3,000 km^ with the atmosphere 
above viewed by 24 fluorescence telescopes placed in 4 buildings. 
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Each WCD consists of a rotomolded polyethylene tank with an area of 10 and 1.55 m high, 
enclosing a Tyvek® liner filled with 12,000 1 of high purity water. The Cherenkov light is detected 
by three large 9" photomultipliers (PMTs) viewing the water tank from the top. Each detector 
is an autonomous unit. The signals from the PMTs are processed by the local electronics. The 
power of the whole system is provided by two batteries connected to two solar panels. The time 
synchronization is done by commercial GPS receivers. There is also a communication system 
between each detector and the Central Station of the Observatory. 

The UHECRs are detected through the particle cascade they produce in the atmosphere. With 
the Auger surface detector (SD) we determine their arrival direction from a fit to the particle 
arrival times on the ground. 

During the development in the atmosphere, the air shower increases in size up to a maximum. 
The bulk of the particles reaching the ground are grouped in a thick (hundreds of nano-seconds 
to a few micro-seconds, e.g. 100 m to 1000 m) nearly spherical pancake. The thickness depends 
on the distance to the shower axis, while the radius of curvature of the sphere that can be used to 
describe the shower front approximately depends on the altitude of the shower maximum. 

The precision in the cosmic ray arrival direction achieved from the SD reconstruction depends 
on the clock precision of the detector and on the fluctuations in the first particle arrival time at 
ground. These fluctuations increase with the shower thickness and decrease with the particle 
density. Therefore, the shower front time will be measured with less precision for a detector 
located further away from the core, where particle densities are small and the shower thickness is 
large, than for the same detector closer to the core of the same shower. 

To estimate the primary direction of the primary cosmic ray with the maximum precision, one 
must, among other things, adequately model the measurement uncertainties of each individual 
tank participating in the event. These measurement uncertainties depends on the shower proper- 
ties at each tank location. With an adequate model for the individual measurements, one is able 
to weight the contribution in the determination of the arrival direction correctly. 

The model of the shower front used in the minimization procedure, be it spherical, parabolic, 
or even planar also influences the uncertainty in the arrival direction determination, but not as 
much as the time measurement precision. We will show in section|5]that a parabolic model for the 
shower front adequately describes the data. 

In section |2] we describe a model for the shower front particle distribution from which we 
derive the uncertainty in the time measurement at each tank. In section |3] we implement this 
model for the Auger WCDs. In section|4]we compare it to measurements obtained using two pairs 
of adjacent detector stations located in the Auger surface array. Finally, in section |5] we further 
validate our model by studying the probability distribution of the geometrical reconstruction 
of Auger events. 



2 A model for the measurement uncertainties 

As an estimator for the shower front arrival time (T^) at each detector location on the ground it is 
customary to use the time of the first particle entering the tank. As discussed in the introduction. 
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assuming one has a reasonable description of the shower front geometry, it is the estimation ac- 
curacy on Ts at each tank that will determine the overall precision with which the Auger SD can 
reconstruct the primary cosmic ray arrival direction. 

At a given point on the groimd, we assume that the particle distribution in the shower front can 
be described as a Poisson process with parameter r starting at time to- Of course, both r and to are 
functions of the shower characteristics at the considered ground location. In reality, the particle 
distribution in the shower front is not a uniform Poisson process. At the beginning of the front 
the particle arrival frequency is larger than towards the end. It is however sufficient to assume 
that the frequency is constant over some interval large enough to have a good fraction of the total 
number of particles reaching the ground but less than the total thickness of the shower front at 
that location. 

The first particle arrival time as recorded by the detector clock is then given by: 

Ts = to + Ti (1) 
where Ti follows an exponential distribution function with decay parameter r. 

Note that, by construction Ts (our shower front arrival time estimator) is always larger than 
the true (unknown) shower front time to- is a biased estimator of to with a bias given by the 
expectation of Ti, E[Ti] = t.Ht is known it would be interesting to construct a new shower front 
estimator as T'^ = Tg — r. However, r can only be estimated from the shower parameters measured 
at the ground and is therefore subject to fluctuations. The new estimator has therefore slightly 
larger fluctuations than TsEl For an schematic view of our notation see figured] 

The variance of Ts {V[Ts]) is given by the sum of the detector clock precision (b^) and of the 
variance of the random variable Ti, V[Ti] = r^. In practice r is unknown and must be estimated 
from the data itself. We reverted to use (see section |3ll as an estimator of r the ratio of the time of 
arrival of the k-th particle (T^) to k: f = T^/k. The variance of Ts then become^ 

where 6^ represents the GPS uncertainty and the resolution of the flash analog-to-digital convert- 
ers (FADCs) of the tanks. 



3 Model implementation 

To estimate equation |2] for each tank participating in the event we want to reconstruct, we must 
define Tk and k from the data measured at the tank. For k we choose to use the total number of 
particles entering the tank (n in the following). 

The shower signal is recorded in the Auger WCD by the FADCs, which give the signal values in 
bins of 25 ns as a function of time. From these FADC traces one can build several time parameters 

^ When the radius of curvature of the shower front is estimated during the geometrical minimization (which happens 
for tank multiplicities larger than 5) the bias is absorbed in the radius of curvature. For low multiplicity events with 
only 3 tanks where there is no constraint to determine the geometry, this bias induces a systematic in the pointing but 
much smaller than the pointing resolution itself. 

^V[Ti] is not strictly = (Tfe/fe)^ but f^{k — l)/(fc + 1) because Tk is measured from the data themselves. 
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such as the time it takes to reach 50% of the total integrated signal in the trace: T50 (see figured]). 
Following our remark in section IH we decided to use 2T50 for T„ rather than Tioo to account for 
the fact that the Poisson process is not uniform over the whole shower front. 

To calculate the total number of particles, n, we assume that all particles hit the detector with 
the same direction as the shower axis, and that the muons (which arrive first) are the ones that 
contribute most to the time measurements [3 J. Then, we obtain n as the ratio between the total 
integrated signal in the FADC traces in units of vertical equivalent muons [4] (the signal deposited 
by a central and vertical downgoing muon in a tank), divided by the average muon track length 
TL{9) of a shower with zenith angle 0, normalize to the vertical height {h): 

SvEM 



The average track length TL{6) can be computed as the ratio of the detector volume {V) and 
the area {A) subtended by the arriving particles at zenith angle 9. For the Auger cylindrical tanks, 
it is given by: 

rT((}\ = ^ = '"'"^ (A\ 

^' A 7rrcos(e) + 2/isin(0)' ^' 
where r = 1.8 m is the detector radius, and h = 1.2 m is the water height in the detector. 

We further introduce a scale parameter a to account for the various approximations made to 
build V[Ts] and finally define it as: 

2/^2 Tso \ ^ n — I 2 



' n J n + 1 



Parameters a and b can be determined from the data of the two pairs of adjacent stations avail- 
able in the Auger array (see section 141]) . We expect, however, a to be close to unity and to b be 
given by our GPS system resolution ||5l (10 ns) and the resolution of our FADC traces (25/ \/l2 ns), 
giving 6 ~ 12 ns. 



4 Comparison with data 
4.1 Fitting the doublet data 

Two pairs of adjacent surface detector stations separated by 11 m ("doublets") have been installed 
in the field of the Auger Observatory. These pairs enable comparison of timing and signal accuracy 
measurements . 

To adjust the constants a and b we used all the T4 [6] events from April 2004 to the end of 
December 2006, with at least one doublet information available. We defined the time difference 
AT = dTi - dT2 where dTi {dT2) is the time residual of the first (second) twin tank to the fitted 
shower front. Note that using the residual difference assures that AT does not depend on the 
global shower front shape, since the twin tanks are only 11m apart. 
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The knowledge of AT only requires the knowledge of 9 and (p with moderate precision (a few 
degrees is enough). We further required | AT | to be smaller than 200 ns to avoid the tails of the 
distribution and also a good overall reconstruction. With these cuts we had 2037 events to fit our 
two parameters a and b. 

Since the distribution of AT/ y^F[AT] is close to a Gaussian (see figure|2l), we fit the parameters 
by maximizing the following likelihood function (L): 

where is the total number of events and yfATfc] = FfTi fc] + y[T2ji;] is calculated from the 
variance of Ts for each member of the doublet and for each event k using equation|5l Maximizing L 
is equivalent to minimizing Z = — 2 ln(L), which was done using the Minuit package. We obtained: 

a2 = 0.97 ±0.04, 
P = (132 ± 11) ns^. 

These results are in remarkable agreement with the values expected from the model definition: 

= 1 and 5^ = 144 ns^. 

As mentioned earlier, we expect the distribution of the variable AT to be essentially insensitive 
to the details of the geometrical reconstruction used to derive it. This is because the two tanks 
in a doublet are very close together and AT can be computed with moderate precision based on 
the shower angles (a few degrees) and regardless of the general structure of the shower front. 
We verified that indeed the parameters a and b did not depend on the nature and quality of the 
initial geometrical fit used when analyzing the doublet data. We used, for example, values for 
a and b different from 1 and 12 during the shower geometrical fit, and also replaced our time 
variance model by a constant, we always obtained the same final parameters values (within the 
uncertainties). This shows that using the doublets and the estimator AT we are really independent 
of the details of the geometrical reconstruction used to determined the model parameters. 



4.2 Model quality 

The value of L at the minimum (10.9 per degree of freedom) is quite meaningless (unlike a 
value), therefore we need to evaluate the quality of the fit using alternative tests. 

We computed the quality = Yl^=i v^f^f' which should be close to the number of de- 
grees of freedom if y[ATfc] properly models the variance of AT. We obtained at the minimum 
X^/ndof = 1.00097. 

We also generated a Monte Carlo sample based on the doublet data. For each event with a 
doublet we used the values of T50 and signal {S) of the doublet members and draw a new time 
value according to our shower front model, plus a 25 ns uniform distribution corresponding to 
our FADC bin size and a Gaussian distribution with o" = 10 ns corresponding to the GPS clock 
accuracy. We then calculated AT as the difference of these two times and reproduced our analysis 
on this MC sample. 
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From the minimization of Z, we obtain on the MC sample: 

= 0.74 ±0.05, 
62 = (215 ± 28) ns^. 

with a value of Z = 10.6 and = 1.0011 at the minimum, which are similar to what we obtained 
with the real data and also corresponds to the injected parameters. 

In addition, in figure|3]we plot for comparison AT distributions for the doublet data (red-solid) 
and for the Monte Carlo (blue-dashed). As can be seen, the agreement is good. 



5 Validation 

If the time variance model describes correctly the measurement uncertainties, the distribution of 
AT/^y [AT], where V[^T] = V[t[^'^] + l^[rf^], should have unit variance. 

We show on figure Hthe RMS of the distributions of /^T/^JV[l\T] for various bins in cos(0) 
(top), in signal (middle), and in distance to core (bottom). In the three cases, it is almost constant 
and close to unity, which shows that our model for the time uncertainty is in good agreement with 
the experimental data. It is important to remark that our model does not depend directly on the 
distance of the tank to the shower core. Therefore, the results obtained in this case can be regarded 
as an additional proof of validity. 

We also studied the distribution of the probability of the T5 events HH with 4 or more sta- 
tions. In figure|5]we show this distribution for all events (top), for events with zenith angle smaller 
than 55° (middle) and events with zenith angle larger than 55° (bottomjl. For the three cases, the 
distributions are almost flat as it should be in the ideal case. It is important that the flatness is ob- 
served both for large and small zenith angles separately, which means that the model for the time 
uncertainty works well for all angles without compensating one set from the other. The flatness 
of these distributions indicates that our model properly reproduces the experimental uncertainty. 
It also indicates that the shower front model used in our fit (parabolic front with adjusted radius) 
is adequate. 

Finally, to illustrate the sensitivity of the probability distribution to the value of the variance 
used in the fit, we plotted in figure |6] the distribution obtained when one artificially varies the 
overall variance (e.g. normalization. In black (solid), the original curve, in red (dashed) for 
a variance reduced by 20% (individual mutliplied by 1.2), in blue (dotted) increased by 20% 
(individual multiplied by 0.8). This 20% in is equivalent to a 10% change in a. We therefore 
conclude that we properly reproduce the measurement uncertainties within 10% level. 

^We only plot probabilities larger than 1% to avoid the large peak at zero corresponding to badly reconstructed 
events, which corresponds to about 9% of the total data. 
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6 Conclusions 



We developed a model to describe the measurement uncertainties of the arrival direction of the 
shower front for an air-shower surface detector array. This model is used for the Auger WCD and 
is based on the shower zenith angle, integrated signal and rise time measured in the tanks. Ab- 
solute predictions of our model parameters for the time uncertainty measurements are consistent 
with the ones obtained from analysis of adjacent tank data. 

We performed numerous tests to validate both our estimation methods and our modeled val- 
ues. All cases indicated a correct procedure and a good agreement with the experimental data. 

This model together with the proper description of the front geometry will allow for an optimal 
determination of the shower arrival direction. Properly taking into account our measurement 
uncertainties also allows for an event by event determination of the precision on the reconstructed 
arrival direction. 
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Figure 1: Example of an integrated FADC trace for a detector at 700 m from the shower core 
position. Tg is the shower front arrival time, Ti is the arrival time of the first particles, and T50 is 
the time it takes to reach 50% of the total integrated signal in the trace. 




Figure 2: Distribution of AT/,/lV[AT]) where V[AT] = V[Ti] + V[T2] with Ti {T2) calculated using 
equation |5] for the first (second) twin tank, with a = 1 and b = 12 ns. 



8 



Entries 2037 
Mean -0.6001 
RMS 68.20 




Time difference (ns) 



Figure 3: Comparison of the AT distributions for doublet data (red-solid) and for MC events (blue- 
dashed). 
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Figure 4: The RMS of the distribution of AT/ y^V[AT], as a function of the shower zenith angle 
(top), the average signal in the doublet detectors (middle), and the distance to the shower core 
(bottom). 
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Figure 5: The probability distribution for all events (top), events with zenith angle smaller 
than 55° (middle), and events with zenith angle larger than 55° (bottom). In the last figure the 
distribution is plotted with two different scales, the same than the others (full line) for comparison 
reasons and a zoom (dashed line) to see the details. 
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Figure 6: Distribution of the probability of x^, in black (solid) the original one, in red (dashed) 1.2 arid 
in blue (dot) 0.8 x^- 
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